Scott Cole
12 March 2016
The phase of an ongoing neural oscillation is often a useful metric for analyzing neural signals. The common way to obtain a phase time series is to band-pass filter the data in the frequency range of interest (e.g. "beta," 13-30Hz), and then take the angle of the complex signal that is output of the Hilbert transform of the band-pass filtered data.
However, to account for the nonsinusoidal nature of hippocampal theta, numerous authors have calculated instantaneous phase using properties of the raw oscillatory waveform (e.g. Siapas et al., 2005).
In this notebook, we demonstrate the method used in our lab to estimate phase based on the raw waveform.
Imports
In [1]:
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("white")
from misshapen.nonshape import wfpha, findpt, bandpass_default, phaseT
necessary functions
In [2]:
def zerox(x, frange, Fs = 1000.):
# Filter in narrow band
xn, _ = bandpass_default(x, frange, Fs, rmv_edge = False)
# Initialize phase array
L = len(xn)
pha = np.empty(L)
pha[:] = np.NAN
# Find zero crosses
def fzerofall(data):
pos = data > 0
return (pos[:-1] & ~pos[1:]).nonzero()[0]
def fzerorise(data):
pos = data < 0
return (pos[:-1] & ~pos[1:]).nonzero()[0]
zr = fzerorise(xn)
zf = fzerofall(xn)
return xn, zr, zf
def _edgeadd_paseries(amp, fosc, Fs, w = 3):
"""
Undo the removal of edge artifacts done by pacpy in order to align
the extrema with their amplitudes
"""
Ntaps = np.floor(w * Fs / fosc[0])
amp2 = np.zeros(len(amp)+2*Ntaps)
amp2[Ntaps:-Ntaps] = amp
return amp2
Load recording
In [3]:
v = np.load('./exampledata.npy') # DBS setting to plot: before DBS
v = v[:6000]
dt = .001
Fs = 1/dt
t = np.arange(0,len(v)*dt,dt)
flo = (13,30) # Beta band (Hz)
# Visualize only a small part of the recording
t_plot = (3,4) # Time range to plot (seconds)
sampsT = np.arange(len(t))
samps = sampsT[np.logical_and(t>=t_plot[0], t<t_plot[1])]
Extract waveform properties
In [4]:
# Find zero crossings
vbeta, zxr, zxf = zerox(v, flo, Fs)
zxrp = zxr[np.logical_and(zxr>=t_plot[0]*Fs, zxr<t_plot[1]*Fs)]
zxfp = zxf[np.logical_and(zxf>=t_plot[0]*Fs, zxf<t_plot[1]*Fs)]
# Find peaks and troughs
peaksW, troughsW = findpt(v, flo, Fs = Fs, )
peaksWp = peaksW[np.logical_and(peaksW>=t_plot[0]*Fs, peaksW<t_plot[1]*Fs)]
troughsWp = troughsW[np.logical_and(troughsW>=t_plot[0]*Fs, troughsW<t_plot[1]*Fs)]
# Estimate waveform phase
bphaW = wfpha(v, peaksW, troughsW)
# Estimate phase time series using traditional Hilbert-based method
bphaH = phaseT(v, flo, Fs)
Visualize waveform phase estimate method
In [5]:
fig = plt.figure(figsize=(15,10))
plt.subplot(4,1,1)
plt.plot(t[samps]-t[samps[0]], v[samps],'k',linewidth=3)
plt.ylabel('Voltage ($\mu V$)',size=20)
#plt.ylim((-300,700))
#plt.yticks([-300,0,700])
plt.tick_params(labelsize=20)
plt.setp( plt.gca().get_xticklabels(), visible=False)
plt.subplot(4,1,2)
plt.plot(t[samps]-t[samps[0]], vbeta[samps], 'k',linewidth=3)
plt.plot(t[zxrp] - t[samps[0]], np.zeros(len(zxrp)), 'go', ms=8,label='Rising zerocross')
plt.plot(t[zxfp] - t[samps[0]], np.zeros(len(zxfp)), 'yo', ms=8,label='Falling zerocross')
lgd = plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),fontsize=20)
plt.ylabel('Beta-filtered\n Voltage ($\mu V$)',size=20)
#plt.ylim((-400,400))
#plt.yticks([-400,0,400])
plt.tick_params(labelsize=20)
plt.setp( plt.gca().get_xticklabels(), visible=False)
plt.subplot(4,1,3)
plt.plot(t[samps]-t[samps[0]], v[samps], 'k-',linewidth=3)
plt.plot(t[peaksWp]-t[samps[0]],v[peaksWp],'mo', ms=8, label='peak')
plt.plot(t[troughsWp]-t[samps[0]],v[troughsWp],'co', ms=8, label='trough')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),fontsize=20)
plt.ylabel('Voltage ($\mu V$)',size=20)
#plt.ylim((-300,700))
#plt.yticks([-300,0,700])
plt.tick_params(labelsize=20)
plt.setp( plt.gca().get_xticklabels(), visible=False)
plt.subplot(4,1,4)
plt.plot(t[samps]-t[samps[0]], bphaH[samps], color='0.5',linewidth=3, label='Bandpass')
plt.plot(t[samps]-t[samps[0]], bphaW[samps], 'k',linewidth=3, label='Broadband')
plt.ylabel('Beta phase',size=20)
plt.xlabel('Time (s)',size=20)
#plt.ylim([-4,4])
#plt.yticks([-np.pi,0,np.pi],['$-\pi$','0','$\pi$'])
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),fontsize=20)
plt.tick_params(labelsize=20)
plt.tight_layout()
Raw voltage (row 1). Zero-crossings are identified in a signal bandpassed in the frequency range of interest (row 2). Peaks and troughs are then identified by finding the timepoints of maximal and minimal voltage between each extrena (row 3). Instantaneous phase is found by linear interpolation over time between peaks (phase = 0) and troughs (phase = -π,π). Note the difference between the traditional method ("bandpass", gray) and the waveform-based method ("broadband", black).